!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>


!     subroutine LIMBES_init fetch all data necessary for the simulation
!     -> fluid properties in SIMUL_path/
!     -> geometry configuration in SIMUL_path/
!     -> initial distribution function f in SIMUL_path/f_init.dat

subroutine LIMBES_init()
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: ni,vi,nv
  integer*8 :: n1,n2,n3,n4,n5,n6,n7
  real*8 :: r1,r2,r3,r4,r5,r6,r7,r8
  real*8 :: T0
  character*2 :: char,char1,char2
  character*50 :: arg
 
  !     pi=3.14
  pi=4.0d+0*datan(one)
  if (iargc().EQ.0) then
     print*,'You forget to indicate the parameter file !!!'
     stop
  endif
  call getarg(1,arg)
  !     fetching simulation parameters
  !  open(10,file='./simul_parameters.in')
  open(10,file=arg)
  read(10,*) ! time steps number 
  read(10,*) SIMUL_step_number
  read(10,*) ! number of results to save 
  read(10,*) SIMUL_save_number
  read(10,*) ! reference temperature
  read(10,*) T0
  read(10,*) ! gravity
  read(10,*) SIMUL_gx,SIMUL_gy
  read(10,*) ! quadrature order 
  read(10,*) SIMUL_Nquad
  read(10,*) ! Hermite expansion order
  read(10,*) LIMBES_NH
  read(10,*) ! dt/tau ratio
  read(10,*) SIMUL_dttau
  read(10,*) ! transport sheme
  read(10,*) LIMBES_transport_scheme
  read(10,*) ! thermal or isothermal simulation
  read(10,*) SIMUL_thermal
  read(10,*) ! data path
  read(10,*) SIMUL_path
  close(10)

  LIMBES_nv=SIMUL_Nquad**two
  nv=LIMBES_nv   

  !     fetching mesh information
  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'mesh.in')
  read(10,*)
  read(10,*) ! cells number
  read(10,*)
  read(10,*) SIMUL_cell_number, SIMUL_dx,SIMUL_nx,SIMUL_ny
  read(10,*) VTK_xmin,VTK_xmax,VTK_ymin,VTK_ymax
  !------allocation LIMBES_var------------------
  call  init_LIMBES_var(SIMUL_nx,SIMUL_ny,SIMUL_cell_number,nv,SIMUL_save_number)
  call check_init_LIMBES_var()
  !--------------------------------------------
  do ni=1,SIMUL_cell_number
     read(10,*) n1,n2,n3,n4,n5,n6,n7
     LIMBES_cell_xi(ni)=n2
     LIMBES_cell_yi(ni)=n3
     LIMBES_cell_ni(n2,n3)=n1
     LIMBES_neighbour_RLUD(ni,1)=n4
     LIMBES_neighbour_RLUD(ni,2)=n5
     LIMBES_neighbour_RLUD(ni,3)=n6
     LIMBES_neighbour_RLUD(ni,4)=n7
  enddo
  close(10)

  !     fetching fluid properties
  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'fluid_properties.in')
  read(10,*) ! massique fluid constant 
  read(10,*) FLUID_R
  read(10,*) ! collision time
  read(10,*) FLUID_tau
  read(10,*) ! collision time expression
  read(10,*) FLUID_tau_expr
  close(10)

  ! reference temperature and time step
  FLUID_RT0=FLUID_R*T0
  SIMUL_dt=SIMUL_dttau*FLUID_tau

  !     initialazing discrete velocities model 
  call LIMBES_dvm()

  !     fetching boundary data
  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'boundary.in')
  read(10,*)
  read(10,*)
  read(10,*) SIMUL_1wall_cell_number,SIMUL_2wall_cell_number
  !------allocation SIMUL_var------------------
  call init_LIMBES_fluid_config(SIMUL_1wall_cell_number,SIMUL_2wall_cell_number)
  call init_LIMBES_var_wall(SIMUL_1wall_cell_number,SIMUL_2wall_cell_number,nv)
  call check_init_LIMBES_fluid_config()
  !--------------------------------------------
  do ni=1,SIMUL_1wall_cell_number
     read(10,*) n1,char,n2,n3,n4,n5,n6,r1,r2,r3,r4
     SIMUL_1wall_boundary_type(ni)=char
     SIMUL_1wall_cell(ni)=n2
     SIMUL_1wall_RLUD(ni,1)=n3
     SIMUL_1wall_RLUD(ni,2)=n4
     SIMUL_1wall_RLUD(ni,3)=n5
     SIMUL_1wall_RLUD(ni,4)=n6
     SIMUL_1wall_alpha(ni,1)=r1         
     SIMUL_1wall_alpha(ni,2)=r2        
     SIMUL_1wall_alpha(ni,3)=r3         
     SIMUL_1wall_alpha(ni,4)=r4        
  enddo
  read(10,*)
  do ni=1,SIMUL_2wall_cell_number
    
     read(10,*) n1,char1,char2,n2,n3,n4,n5,n6,r1,r2,r3,r4,r5,r6,r7,r8
     SIMUL_2wall_cell(ni)=n2
     SIMUL_2wall_boundary_type(ni,1)=char1
     SIMUL_2wall_boundary_type(ni,2)=char2
     SIMUL_2wall_RLUD(ni,1)=n3
     SIMUL_2wall_RLUD(ni,2)=n4
     SIMUL_2wall_RLUD(ni,3)=n5
     SIMUL_2wall_RLUD(ni,4)=n6
     SIMUL_2wall_alpha1(ni,1)=r1         
     SIMUL_2wall_alpha1(ni,2)=r2        
     SIMUL_2wall_alpha1(ni,3)=r3         
     SIMUL_2wall_alpha1(ni,4)=r4  
     SIMUL_2wall_alpha2(ni,1)=r5         
     SIMUL_2wall_alpha2(ni,2)=r6        
     SIMUL_2wall_alpha2(ni,3)=r7         
     SIMUL_2wall_alpha2(ni,4)=r8   
  enddo
  close(10)

  !     fetching initial distribution function f
  open(10,file=SIMUL_path(1:strlen(SIMUL_path))//'f_init.dat',form='unformatted')
  do ni=1,SIMUL_cell_number
     read(10) (LIMBES_f(ni,vi),vi=1,nv)       
  enddo
  close(10)

  !     few info
  write(*,*) '----------LIMBES running----------'
  write(*,*) 'Domain size',SIMUL_nx,'x',SIMUL_ny
  write(*,*) 'delta t:', SIMUL_dt,'   delta x:', SIMUL_dx
  write(*,*) 'Discrete velocity model:', LIMBES_nv
  write(*,*) 'Hermite order of equilibrium:', LIMBES_NH
  write(*,*) 'Transport scheme: ', LIMBES_transport_scheme
  write(*,*) 'CFL=',LIMBES_vx(1)/(SIMUL_dx/SIMUL_dt)

end subroutine LIMBES_init

